Second generation wave-function thermostat for ab-initio molecular dynamics 
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A rigorous two-thermostat formulation for ab-initio molecular dynamics using the fictitious La- 
grangian approach is presented. It integrates the concepts of mass renormalization and temperature 
control for the wave functions. The new thermostat adapts to the instantaneous kinetic energy of 
the nuclei and thus minimizes its influence on the dynamics. Deviations from the canonical ensem- 
ble, which are possible in the previous two-thermostat formulation, are avoided. The method uses 
a model for the effective mass of the wave functions, which is open to systematic improvement. 
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I. INTRODUCTION 

Ab-initio molecular dynamics H , also called the Car- 
Parrinello method, allows to study the dynamical and 
finite temperature behavior of molecules, surfaces and 
solids with forces from highly accurate first principles 
density functional simulations j| ||. This approach has 
revolutionized the way electronic structure calculations 
are done to date. Two main variants of this approach 
have established themselves. The original approach of 
Car and Parrinello|l| uses a fictitious Lagrangian to de- 
duce a dynamical equation of motion for the wave func- 
tions, while the so-called "exact" Born Oppcnhcimcr 
dynamics performs self-consistency loops for each set 
of atomic positions. 

The underlying idea of the fictitious Lagrangian ap- 
proach, which treats the electronic wave functions as 
dynamical fields that obey a Newton-type equation of 
motion. If the temperature attributed to the motion of 
the wave functions is sufficiently low, the wave functions 
propagate close to the Born-Oppenheimer surface (i.e. 
the instantaneous electronic ground state). Thus the 
electronic wave functions can be propagated with rel- 
atively minor computational effort while the nuclei are 
moving. Self-consistency loops for each step of the tra- 
jectory are avoided. 

A difficulty arises from the requirement that the 
wave functions remain sufficiently close to the Born- 
Oppenheimer surface so that physically correct forces on 
the nuclei are produced. This implies that the wave func- 
tions must remain "cold" , while at the same time the 
nuclei are at the, relatively high, physical temperature. 
Hence ab-initio molecular dynamics simulations are in 
principle non-equilibrium simulations. 

In most cases, the heat transfer between the electronic 
and the atomic subsystems is sufficiently slow as a re- 
sult of the separation between electronic and nuclear fre- 
quency spectrap]. For long simulations, however, one 
has to bear in mind that the nuclear and electronic vari- 
ables will eventually, even though very slowly, approach 
thermal equilibrium, where the wave functions deviate 
from the Born-Oppenheimer surface and the forces act- 
ing on the atoms are unphysical. In order to rigorously 



maintain a stable simulation, two thermostats jfj, [7| are 
introduced M. One thermostat keeps the nuclei at their 
physical temperature and the other one adsorbs the addi- 
tional heat transfer to the wave functions. The optimum 
temperature for the wave-function dynamics, which is re- 
quired so that the wave functions can follow the nuclei 
adiabatically, has been discussed earlier. g] 

However, the previous two-thermostat formulation has 
three deficiencies: 

• For systems consisting of parts that are only in 
weak thermal contact, the thermostats do not guar- 
antee that the canonical ensemble is correctly sam- 
pled. The heat transfer from the nuclei to the wave 
functions is larger for some atoms than for others, 
which may result in different effective temperatures 
for different parts of the system. 

• The target kinetic energy for the wave function dy- 
namics is constant over time and is determined ac- 
cording to the target temperature of the atoms. For 
small systems the fluctuations in the atomic kinetic 
energy can be sizeable. In that case the wave func- 
tion thermostat will heat the wave functions when 
the atomic kinetic energy is low and cool more than 
necessary if the atoms have a high kinetic energy. 
This has an adverse effect on the nuclear trajecto- 
ries. 

• The goal of the wave function thermostat is to keep 
the wave functions cold. The thermostat however 
induces fluctuations in the wave function kinetic 
energy and thus introduces undesired heating se- 
quences. 

In this paper, I analyze the trajectories of the nuclei 
and propose an new formulation of the two-thermostat 
method, which overcomes these deficiencies. The new 
method links the target temperature for the electronic 
wave functions directly to the instantaneous motion of 
the atoms. Thus the wave function thermostat adapts to 
the temperature fluctuations of the nuclei, and minimizes 
its influence on the atomic trajectories. Secondly, the 
indirect influence of the wave-function thermostat on the 
atomic trajectories is compensated using an additional, 



opposing friction term in the equations of motion for the 
atoms. 

A detailed numerical analysis of the errors in the 
forces resulting from the deviations from the Born- 
Oppenheimer surface and due to the previous two- 
thermostat formulation been performed H. 

In Sec. II the previous two-thermostat formulation is 
restated, which allows me to introduce my notation. In 
Sec. Ill effective equations of motion for the nuclei un- 



der the influence of the wave function dynamics are de- 
rived. Sec. IV discusses a potential thcrmodynamical in- 
stability in the previous two thermostat formulation of 
Car-Parrinello dynamics. Sec. [V] describes the new two 
thermostat formulation. An approximate scheme for de- 
riving the parameters of the new theory and technicalities 
of the present implementation are giv en in Sec. N% Test 



calculations are presented in Sec. VII 



II. 



EQUATIONS OF MOTION WITH TWO 
THERMOSTATS 



Let me first describe my notation. The electrons are 
described by one-particle wave functions \^ n )- The wave 
functions |W n ) are related by a linear transformation 
\^n) — T\4f n ) to the variational parameters \^ n )- Of- 
ten the variational parameters are vectors and the trans- 
formation is defined by a basis set \\k) in the form 
l^n) = J2k \Xk)^k,n- Here, I have in mind the Projector- 
Augmented Wave (PAW) method jlQ], where the varia- 
tional parameters are themselves fields, namely the so- 
called pseudo wave functions \^f). The pseudo wave func- 
tions are expanded in a basis set, a detail that does not 
concern here. The pseudo wave functions of the PAW 
method are conceptionally identical to the wave func- 
tions of the pseudopotential formalism. 

In the ab-initio molecular dynamics method with the 
previous two-thermostat formulation, the following cou- 
pled system of equations of motion is solved. 

m*|*„) = -#!#„) +^0|# m )A m ,„-m*|i' n )x* 

m 

MiRi = Fi-MiRiXR. 



Q R x R = 2{Y,\MiR 2 i-\gk B T). 



(v&nl^m) — {^n\0\^ m )- The Lagrange parameters intro- 
duced to keep the wave functions orthogonal are denoted 
as A mj „. The mass of the wave functions is m*. In prac- 
tice an operator diagonal in a plane wave representation 
is used. The nuclear masses are denoted as Mi, where 
i is the index of the corresponding atom. Qq, and Qe 
are the "masses" of the thermostat variables xq, and xr 
for wave functions and nuclei respectively. These masses 
determine the reponse time and the dominant frequency 
of the thermostats, g is the number of nuclear degrees 
of freedom and kg is the Boltzmann constant. Ekinfi is 
a parameter that determines the target kinetic energy of 
the wave functions in the simulations. Its value is cho- 
sen according to an estimate of the Born-Oppenheimer 
kinetic energy of the wave functions. In the entire paper 
Hartree atomic units {% = e = m e = Aneo = 1) are used. 
The present formulation of the thermostat differs from 
Hoover's: the thermostat variable used by Hoover corre- 
sponds to the time derivatives of the thermostat variables 
xr and xq, used here. This choice has the advantage that 
the thermostats obey second order differential equations 
just as the nuclear position and the wave functions. 

Even though the equations of motion given above can- 
not be derived from a Lagrangian formalism, they have 
a conserved energy 

E c = ^(fnlm^l^+^-M^ + EO^),^] 

n i 

+ 7jQ*i| + 2E kinfi xq, + -Qrx 2 r + gk B Tx R .(2) 

The first two terms are the kinetic energies of wave func- 
tions and nuclei, the third is the potential energy, i.e. the 
density-functional total energy. The remaining terms are 
kinetic and potential energies of the two thermostats for 
wave functions and nuclei. 

Simulations with the two thermostats reach a station- 
ary state in which the wave function thermostat vari- 
able xq, exhibits an approximately constant drift to larger 
values, which freezes out the deviations from the Born- 
Oppenheimer approximation. As a consequence of en- 
ergy conservation jll| the atom thermostat variable de- 
creases at an average rate 
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which restores the energy absorbed by the wave-function 
thermostat. 



If E[R, |W)] is the Kohn-Sham total energy functional, 
the pseudo Hamiltonian H is defined such that d P , — 

H\^> n ) and the forces Fi are the partial derivatives of 
the Kohn-Sham total energy functional with respect 
to the atomic positions. The pseudo overlap operator 
O = T^T is obtained from the transformation between 
the true wave functions |\l/) and the pseudo wave func- 
tions \*&). The overlap between two wave functions is 



III. EFFECTIVE EQUATIONS OF MOTION 
FOR THE NUCLEI 

In order to understand the motion of the nuclei with 
the wave functions and thermostats tied to them, let me 
derive effective equations of motions for the nuclei. It is 
useful to consider the atoms as quasiparticles consisting 
each of a nucleus and the wave function cloud, i.e. the 



distortion of the surrounding electron gas. Just as an 
electron distorts a surrounding crystal lattice to form a 
polaron, here a nucleus distorts the electron gas to form 
a quasiparticle called an atom. The distortion of the 
electron gas, respectively its wave functions, will be called 
wave function cloud. The effective equations of motion 
are obtained from the equations of motion in Eq. PP by 
constraining the wave functions to be identical to the 
exact Born-Oppenheimer wave function. The details of 
the derivation are given in App. [M Here only the result 
is shown: 



of a subsystem A imposed by the thermostats is accord- 
ing to Eq. § 



E A = — 2, MiRiin — 2_j Ki.jRiRji^, 



(6) 



Assuming that subsystem A is in thermal equilibrium 
at a temperature T A , the thermal average of RiRj is 
(RiRj)t = Si.jksTA/Mi. Furthermore, the drifts of the 
two thermostat variables are related by Eq. ||.[|l3| Thus 
Eq. H can be transformed to 



Y^MiSu + KiJRj = F t - M t R t i R ~ J2 KijRji* (E a )t = gAk B T A (i*){^^- - — £ ^] , (7) 
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Ki.j is the effective mass tensor as function of the nuclear 



positions 



*U = 2 £<-fl5H"»»|-S57>. (5) 



where \^n°) are the Born-Oppenheimer wave functions. 
Let me discuss this equation: 

• The nuclei obtain an effective mass MiSij + Kij: 
The atoms appear to be heavier than the nuclei 
alone, because also the wave functions need to be 
accelerated, whenever the nuclei accelerate. This 
effect has been realized before [[10| [t2| and correc- 
tions are in common use today. 

• The third term on the right-hand side of Eq. |J is 
related to the friction imposed by the wave func- 
tion thermostat on the wave functions. As the 
thermostat removes kinetic energy from the wave 
functions, also the atoms are cooled down. In par- 
ticular for two loosely coupled subsystems, where 
one subsystem has a different average effective wave 
function mass than the other, a drift of the ther- 
mostat variables cools one subsystem at the cost 
of the other. This may cause deviations from the 
canonical ensemble as discussed in more detail in 
the following section. 

• The last term in Eq. describes the effect of the 
changes of the effective mass as the atoms are mov- 
ing around. This term will not be discussed in more 
detail, because in this paper an approximation with 
the effective masses are independent of the atomic 
positions will be employed, where this term van- 
ishes. 



IV. DEVIATIONS FROM THE CANONICAL 
ENSEMBLE 

Here, the deviations from the canonical ensemble men- 
tioned above are investigated. The rate of energy change 



^ ttA Mi 



where g and gA are the number of degrees of freedom of 
the total system and subsystem A. 

If the thermal average of the Born-Oppenheimer wave 
function kinetic energy, 
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is chosen for Eki n ,o> one obtains 
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where g B — g — g A is the number of degrees of freedom 
in the subsystem B which together with subsystem A 
makes up the complete system. 

Hence, the effect of the thermostats can be described 
as a heat flow from one system to the other, with a heat 
transport coefficient 
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Note that \ x is negative. Thermal equilibrium is only 
reached when total heat transport coefficient among sub- 
systems is positive. This clearly shows that there is a 
thcrmodynamical instability when the physical thermal 
coupling is smaller (in absolute values) than the implicit 
coupling via the thermostats. 

The conditions for this thermodynamic instability are 
that 

• the drift in the wave function thermostat is suffi- 
ciently rapid, 

• Ki^/Mi is very different for the atoms of one sub- 
system as compared to the other, and 

• the thermal coupling between the systems is suffi- 
ciently small, so that the effect of the thermostats is 
appreciable compared to the physical thermal cou- 
pling. 



While these conditions rarely corroborate, already the 
fluctuations, which result from a nearby instability, can 
render the results of a simulations inaccurate. 



V. EQUATIONS OF MOTION WITH 
CONTROLLED HEAT-BACK-FEEDING 

To cure the problems of the two-thermostat formula- 
tion, terms are added to the equations of motion that 
compensate the effect of the finite effective mass of the 
wave functions in the down- folded equations of motion for 
the atoms. I proceed with the assumption that the effec- 
tive mass tensor of the wave functions is known. Later, 
I will discuss an approximate expression to be used in 
practice. 

I propose a new set of equations of motion, which is 
the main result of this paper: 

m 

^2(MiS itj - Ki^Rj =Fi- MiRiXR + Y^ KijRjXv 
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Q R x R = 2{Y J \M l R 2 l - l -gk B T). (11) 



The step function 9(x) in the equation of motion for the 
wave function thermostat will be discussed in detail be- 
low. Even if the step function is removed, these equa- 
tions create a stable and energy conserving dynamics. 
Its role is to shut the thermostat down for those periods 
of time, where the thermostat would otherwise heat the 
wave functions. 

The system of equations Eq. [ll] has a conserved energy 



n i,j 

-Qyx\ + -Qrx\ + gk B Tx R 



(12) 



The second line of Eq. U2 is the kinetic energy of the 
free motion of the wave functions representing deviations 
from the Born-Oppenheimer surface. 

The new system of equations of motion differs in four 
points from the ones used previously. 

• The nuclear mass is reduced by the effective wave 
function mass, so that the effective mass of atoms, 
namely the sum of the reduced mass of the nuclei 
and mass of the wave function cloud add up to the 



true nuclear mass. This correction has been previ- 
ously suggested [|10| |l^] and is common practice in 
state-of-the-art ab-initio molecular dynamics simu- 
lations. 

• A friction force ^ ■ KijRjiy acting on the atoms 
has been added. This term is controlled by the wave 
function thermostat and opposes the drag from the 
wave functions, which themselves are slowed down 
by the wave function thermostat. 

The rate of energy added to or removed from the 
physical system by the wave functions thermostat 

is E„(*«l TO y>l*n) -Y,i,j K i,J R i R j}%*- This en- 
ergy transfer is proportional to the deviation of the 
wave functions from the Born Oppenheimer sur- 
face. Hence the thermostat acts only on the free 
oscillations of the wave functions. 

• The estimate | ^ KijRiRj for the instanta- 
neous Born-Oppenheimer kinetic energy has been 
introduced into the feedback equations for the wave 
function thermostat instead of a constant target en- 
ergy used previously. This modification was neces- 
sary in order to obtain a conserved energy. 

Besides restoring energy conservation, this term 
has another beneficial role, which is most important 
for systems with large kinetic energy fluctuations. 
The wave function thermostat adapts to the instan- 
taneous kinetic energy of the nuclei and acts only on 
the deviations of the wave functions from the Born- 
Oppenheimer surface: The wave function thermo- 
stat would remain inactive for Born-Oppenheimer 
wave functions. 

• The wave function thermostat shuts down in an 
energy conserving fashion when it would otherwise 
heat the wave functions. This is discussed in the 
following. 

Despite the more complex equations of motion the con- 
served energy has one term less. There is no potential 
energy term related to the wave functions, which has im- 
portant implications. The thermostat variable oscillates 
repeatedly between cooling and heating the wave func- 
tions. Heating the wave functions is not desirable, as 
one would like to keep them at their lowest temperature 
compatible with the dynamics. 

A better strategy is to couple the wave function ther- 
mostat to the system only when cooling is required. This 
can be done in an energy conserving way, since the con- 
tribution of the wave function thermostat to the con- 
served energy vanishes whenever x^ vanishes. Thus the 
equation of motion for the wave function thermostat in 
Eq. O has been modified by introducing a step func- 
tion 9(xqi). 6(x) is the Heaviside step function defined as 
0{x > 0) = 1, 0(0) = \ and 6{x < 0) = 0. The discretiza- 
tion of the equation of motion for the wave function ther- 
mostat including the step function is not straightforward. 



This equation means the following: The thermostat dy- 
namics is switched on with zero velocity, when the wave 
function kinetic energy grows larger than the target, x^ 
grows and transfers energy from the wave function dy- 
namics into the nuclear subsystem until the velocity of 
the thermostat variable vanishes again. At this point 
the wave function kinetic energy is below its target. In- 
stead of allowing the thermostat to heat the system, the 
thermostat variable x^, is instead kept constant, imply- 
ing that the thermostat is effectively switched off. If the 
equation of motion is solved continuously, the velocity 
of the thermostat vanishes exactly, when the thermo- 
stat shuts down. The velocity will remain zero, until 
the wave function kinetic energy grows above its target, 
and thus does not affect the dynamics for that period of 
time. The thermostat switches on again, when the ki- 
netic energy grows above the target. Thus the trajecto- 
ries proceed without any perturbation as long the wave 
function kinetic energy remains below its target value. 
In a discretized equation of motion the velocity of the 
thermostat must be explicitely reset to zero, whenever it 
would otherwise turn negative. 

To summarize, the thermostat resets the kinetic en- 
ergy to a lower value if the wave functions become too 
hot, while transfering the energy into the nuclear dynam- 
ics. Note that even when the wave function thermostat 
is on, the atomic trajectories are affected only if the ef- 
fective mass tensor is inaccurate. Hence the thermostat 
can operate more strongly before affecting the atomic 
dynamics. 

The thermostat variable of the atom thermostat does 
not experience any longer a steady drift in the new formu- 
lation. This can be deduced from the conserved energy 
expression as follows: A drift of the thermostat variables 
is possible in principle, because a fixed translation of a 
thermostat variable does not affect the dynamics and is 
therefore not observable. The time averaged drift of Xr 
vanishes, because the potential energy of the atom ther- 
mostat is the only term in the conserved energy that de- 
pends on the thermostat variable. All other terms in the 
conserved energy are observable and therefore stationary. 

A drift of the wave-function thermostat variable does 
not affect the conserved energy. Hence, a drift of x^ is 
possible. This drift counter-balances the heat flow from 
the nuclei to the wave functions but does not add to the 
conserved energy. 

In the remainder of this section I will discuss a pos- 
sible dynamical instability and how it can be avoided. 
The Born-Oppcnhcimcr kinetic energy enters with a neg- 
ative sign, which may point to a possible instability of the 
wave functions, when the bare mass-tensor with elements 
MiSij —Kij turns out not to be positive definite: A par- 
ticle with mass accelerates in the opposite direction of the 
forces, and runs away from a minimum. The condition 
for the bare mass tensor to be positive definite poses a 
strict upper bound on the wave function mass m^. 

In practice one needs to be even more restrictive: 
There is an internal excitation of the quasiparticle 



"atom" , where the nucleus oscillates with high frequency 
about the wave function cloud. This mode can be identi- 
fied clearly when the wave functions are frozen and only 
the bare nuclei are allowed to move. If the bare mass of 
the nuclei is small, these oscillations may have high fre- 
quency that require a small time step in the discretiza- 
tion. While this mode falls into the class of deviations 
from the Born-Oppenheimer surface, the thermostat does 
not cure the problem since the wave function kinetic en- 
ergy is smaller than its target | V\ KijRiRj- To avoid 
that problem I suggest to keep the wave function mass 
m^ sufficiently small, so that Ki : i < hMi. 

VI. IMPLEMENTATION 
A. The model of an infinitely dilute gas 

The new equations of motion require an analytic ex- 
pression that approximates the effective mass tensor. As 
a start, I adopt the model of an infinitely dilute gas of 
atoms, which has been used in the context of the pre- 
vious wave function thermostat and mass renormaliza- 
tion. The model of an infinitely dilute gas assumes that 
the Born-Oppenheimer wave functions can be divided up 
into purely atomic contributions, and that those are iden- 
tical to the wave functions of the corresponding isolated 
atoms. 

The rationale for the model of an infinitely dilute gas is 
that the most rapid variations of the wave function occur 
near the nucleus, and are little affected by the bonding 
environment. Thus, already the isolated atom will cap- 
ture most of the relevant contributions. One can envisage 
better approximations for the effective mass of the wave 
functions, which depend explicitely on the atomic posi- 
tions, and which may be derived from a tight-binding like 
description. 

I insert the atomic wave function into the expression 
Eq. for the effective mass tensor of a given atom and 
obtain 



Kij = -2^(# n |Vim*V 3 #„) 



(13) 



It can be shown that in this approximation K% t j is 
diagonal for each atom with three identical matrix el- 
ements on its main diagonal: For symmetry reasons the 
acceleration of an isotropic atom is parallel to the force, 
that is F — XR with some parameter A. Hence, the 
acceleration is an eigenvector of the effective mass ten- 
sor (Mdij + Kij)Rj — XRj. Since A is independent of 
the direction of the applied force, K has three identical 
eigenvalues. Therefore Kij has the form of a unity ma- 
trix times a constant and the identity K^j — 5ij iTrK 
holds. 

Thus one can obtain a more simple form of our model 
effective mass as 



2 



(14) 



In the special case of a G-independent massjlj], the 
weight of the wave function cloud is directly related to the 
kinetic energy of the pseudo wave functions, as reported 
previously Q. 
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A G-dependent wave function mass, usually a mass- 
tensor diagonal in reciprocal space with elements depend- 
ing on the reciprocal space vector of the augmented plane 
waves, is nowadays common practice. It allows to con- 
trol the rapid oscillations of plane waves with large wave- 
and avoids instabilities that occur otherwise 
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vectors 

when the basis set is increased. There are several choices 
for the G-dependence of the wave function mass. I use 
an expression for the effective mass 



m*(G,G') = m%(l + cG 2 )5, 



G,G' 



(16) 



In order to obtain the effective mass tensor, I determine 
the pseudo wave functions of the atom and transform 
them into G-space via a Bessel transform. Then I evalu- 
ate 



A = £/„ fdGG 2 \* n (G)\ 
B = Y,fnfdGG 4 \^ n (G)\ 



(17) 
(18) 



which are combined with the chosen parameters for the 
wave function mass to the effective mass tensor of the 
wave functions 



Ki.4 = §. 



i,j rWf (A; 



cBi). 



(19) 



The variables /„ are the occupation numbers of the one- 
particle states. Typical values for Ai and Bi are given 
in Tab. 0. Note, that these values depend on the choice 
of pseudo wave functions and are not transferable. They 
are listed here solely to provide the reader with a feeling 
of the order of magnitudes involved. 

The effective wave function mass must be substantially 
smaller than the nuclear masses, i.e. ^m^(A + cB) < M, 
in order to obtain a stable dynamics. If this requirement 
is violated the reduced mass of the atom is negative, and 
atoms accelerate in the opposite direction of the force 
acting on them. 



B. Discretized equations of motion 



Atom 


A 


B 


10 4 ^- 

±U M 


10 4 — 

1U M 


H 


0.830 


2.830 


4.518 


15.403 


He 


3.505 


17.128 


4.803 


23.475 


C 


5.837 


22.292 


2.666 


10.181 


o 


14.569 


86.109 


4.995 


29.525 


F 


20.304 


143.785 


5.863 


41.518 


Si 


3.327 


6.320 


0.650 


1.234 


CI 


12.821 


41.841 


1.984 


6.474 


Fe 


10.950 


71.366 


1.075 


7.010 


Ru 


17.606 


95.952 


0.956 


5.208 


Os 


16.874 


86.002 


0.848 


4.320 



TABLE I: Coefficients A and B for Eq.|l9| for different ele- 
ments and their values relative to the nuclear masses. These 
values are not transferable. 



canonical ensemble by tuning the friction via f — x with 
time. 

The equations of motion discretized with a time step A 
are obtained by replacing derivatives by the differential 

quotients x = (x{t + A) - x(t - A))/(2A) + <3(A 2 ) and 
x = (x(t + A) - 2x(t) + x{t - A))/A 2 + 0(A 2 ) as 



x(t + A) 
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-At) - 
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i 



—x(t - A) + F(t) — 



1 



1 + a 

(21) 

where a = /A/2. 

The choice a — yields energy conserving trajectories 
and a = 1 results in steepest descent dynamics. Interme- 
diate values of a are the regime of friction dynamics. 

Because the thermostat can be propagated only with 
the knowledge of the instantaneous kinetic energy, which 
in turn depends on the propagated value of the thermo- 
stat, I extrapolate the thermostat variable for t + A from 
the present and the previous two thermostat values 

x(t + A) = Ax(t) - 6x(t - A) 

+ Ax(t - 2A) - x{t - 3A) + 0(A 4 ) . (22) 

This expression introduces errors of the forces of order 
A 2 , which is consistent with the overall accuracy of the 
Verlet algorithm. 

The step function for the thermostat is implemented 
by resetting x^ to zero for the propagation of wave func- 
tion and nuclear positions, whenever it has a negative 
value otherwise, and by setting x^{t — nA) = x(t), when- 
ever the velocity x^ (t) would become negative during the 
propagation of the thermostat variable. 



All equations of motion are implemented using the Ver- 
let algorithm. An equation of motion for a general coor- 
dinate x has the form 



mx = F — mxf 



(20) 



where / is a friction coefficient, which may be constant 
or imposed by a thermostat. The thermostat x creates a 



VII. TESTS 

In order to test stability and accuracy of the method I 

investigated two systems. A simulation of carbon monox- 
ide shall illustrate how the thermostats adapt to large 
fluctuations of the nuclear kinetic energy. This test case 
allows also a test of the accuracy. Iron has been used 



as example for a metal and shows the stability of the 
trajectories against frequent band crossings. 

The simulations described in the following have been 
performed with the Perdew-Burke-Ernzcrhof density 
functional flq, 17 1. The plane wave cutoff for the wave 



functions has been set at 30 Ry and the one for the den- 
sity at 60 Ry. A time step of A = 10 a.u.= 0.12 fs has 
been used together with a G-independent wave function 
mass of 1000 m e . The frequency of the thermostat for 
the nuclei has been chosen to 10 THz, which is unusually 
large. The frequency of the wave function thermostat has 
been set to 100 THz. 



A. Carbon Monoxide 

The unit cell size is a 13 A fee unit cell. The elec- 
trostatic interaction between periodic images has been 
subtracted. Q A lxs+lxp+lxd projector set has been 
used, that is one projector for every relevant set of angu- 
lar momentum quantum numbers. The one-center den- 
sity has been expanded into spherical harmonics up to 
angular momentum of I = 2. Parameter A = 5.87 for 
carbon and A — 15.42 for oxygen have been used, re- 
sulting in an effective mass of the wave function cloud 
amounts to 17 % and 35 % of the physical nuclear mass 
for carbon and oxygen respectively. Rotations and trans- 
lations have been frozen out by constraints so that a tru- 
ely one-dimensional system is studied. The simulation 
has been performed at 1000 K. 

The experimental value for the CO stretch vibra- 
tion is 2170 cm -1 .p9[.The experimental bond-length is 
1.1283 A. [P0|1 Fixed point calculation predict a bond- 
length of 1.138 A and a stretch frequency of 2125 cm -1 . 
These deviations are in the range of errors expected for 
the density functional used. 

Simulations lasting several picoseconds have been 
performed. pl[ A sequence of 0.3 ps is shown in Fig. [jL 
The deviation of the wave-function kinetic energy from 
the Born-Oppenheimer surface is -0.05 times the vari- 
ation of the potential energy and overall smaller than 
5 meV. This indicates that the dilute atomic gas model 
overestimates the effective masses by about 10%. On the 
other hand 90 % of the error has been removed. 

The reduced mass of atoms with the wave function 
cloud is 25 % of the true reduced mass, which, without 
correction, would result in an overestimate of 10% in the 
frequencies. Given the error in the effective masses one 
can expect frequencies that overestimate the frequency 
derived from static calculation by 1 %, which is good 
agreement with the 2148 cm' 1 obtained from averaging 
vibrational periods during five picoseconds. 

The friction measured in units of %■ remains below 

4 x 10~ 3 -j-, which is very small indicating that the heat 
transfer from wave functions to the nuclei is small in 
this system. The total energy is conserved to within 

5 meV/ps and scales down with the size of the time step. 
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FIG. 1: Energies of carbon monoxide versus time in ps. Top: 
non-Born-Oppenheimer kinetic energy of the wave functions 
in meV (full line) and total energy in units of O.leV (dashed 
line) displaced vertically. Middle: instantaneous "temper- 
ature" in Kelvin. Bottom: friction imposed by the ther- 
mostats in units of 2/A for wave functions (full line) and 
nuclei (dashed line). 



B. Iron 

As a test of a metallic system I have chosen 7-iron. 
An 8 atom fee unit cell has been used. Since only the 
T-point was included in the k-point sampling, I do not 
expect the simulation to be a realistic description of the 
material. The simulation temperature is 1185 K, at the 
martensitic phase transformation temperature. The pa- 
rameter A was about 12 % larger than the kinetic energy 
per atom, in order to account for the promotion of s- to 
d-electrons as one goes from an atom to the solid. The 
effective mass of the wave function per atom is 4.75 a.m.u 
(a.m.u=M( 12 C)/12), about 9 % of the nuclear mass. The 
atom thermostat had a period of 0.1 ps, and the wave 
function thermostat had a period of 0.01 fs. The band 
gap due to finite k-point sampling is typically about 0.1- 
0.2 eV. 

The results for the simulation of one picosecond is 
shown in Fig. 0. The total energy variation is of or- 
der 1 eV. The mean average kinetic energy related to the 
non-Born Oppenheimer motion is 3 meV. The total en- 
ergy drifts with 0.47 meV/ps. The typical deviation from 
the Born-Oppenheimer surface is 10-15 meV. 

Most of the time the thermostat is switched off as the 
kinetic energy remains below the target defined by the 
effective masses. The most pronounced quenching se- 
quence occurs at 0.7 ps. Here the sharp increase of the 
wave function kinetic energy indicates a band crossing. 
A band crossing results in a randomization of the wave 
functions as the occupied state changes its character into 
that of the formerly unoccupied state, and would other- 
wise render the remainder of the simulation useless. The 
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FIG. 2: Energies of austenite (fee-Iron) in an 8 atom super- 
cell versus time in ps. Top: total non-Born-Oppenheimer 
kinetic energy of the wave functions in meV (full line) and to- 
tal energy in units of 0.1 eV (dashed). Middle: instantaneous 
"temperature" in Kelvin. Bottom: friction imposed by the 
thermostats in units of 2/A for wave functions (full line) and 
nuclei (dashed line). 



wave function thermostat brings the wave functions back 
to the Born Oppenheimer surface, while the perturbation 
of the nuclear dynamics during this rather strong quench 
is minimized by the opposing force acting on the nuclei. 



VIII. CONCLUSION 

In summary, a new formulation of the two thermostat 
approach for ab-initio molecular dynamics has been pre- 
sented. The approach aims at controlling only the devi- 
ations from the Born-Oppenheimer wave functions. The 
influence on the Born-Oppenheimer motion of the wave 
functions and the nuclear motion is minimized by addi- 
tional forces opposing the indirect friction of the atoms. 
Furthermore the thermostat is active only if the wave 
function kinetic energy grows beyond its estimated Born- 
Oppenheimer value. The new two-thermostat formula- 
tion can be applied to small systems with large fluctua- 
tions of the nuclear kinetic energy and the fictitious Born- 
Oppenheimer wave function kinetic energy. 

The new approach rests on an expression for the effec- 
tive mass tensor of the wave functions. A simple formula 
has been derived from the previously employed model of 
an infinitely dilute gas. Systematic improvements of ef- 
fective mass tensor, which will improve the quality of the 
simulation, can be envisaged. 
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APPENDIX A: DOWNFOLDING THE WAVE 
FUNCTION DYNAMICS 

Here, the effective equations of motion Eq. ^ for the 
nuclei are derived, which include the forces of the wave 
functions acting on the nuclei. 

The starting point is the following set of equations 



m*|\T/„) 
MRi 



-H\V r . 



6\^ m )-m^ n )x^ 



(Al) 



The wave functions are constrained to remain exactly 
on the Born-Oppenheimer surface, 



tf n (r,t) = *£"(r,i2(t)) 



(A2) 



The Born-Oppenheimer wave functions ^>^°(r, R4) are 
the ground state wave function for a given set of atomic 
positions Ri. 

Because the Born-Oppenheimer wave functions depend 
on the nuclear positions, the forces acting on the wave 
functions translate into additional forces acting on the 
nuclei. Thus effective equations of motion for the atoms 
are obtained. The atoms are now "quasiparticles" con- 
sisting of nuclei and the wave function clouds following 
them. 

The constraints are enforced by the method of La- 
grange multipliers: The constraint forces, which describe 
the effect of the wave function cloud, are the derivatives 
of a "constraint energy" 
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with the auxiliary fields <!>„ acting as Lagrange parame- 
ters. The resulting constraint forces are 



F; 



l*„ 



F 



(*n| 



Ri 



dE c 

dE" 
d|*„> 



dE c 
" dBi 



\®n) 



9* 



BO 



dR,, 



|$n) + ($, 



a* 



BO 



dR,, 



(A4) 
(A5) 

-£A6) 



The constraint forces are inserted in the equation of 
motion for the wave functions 






m 

(A7) 
With the help of the constraint condition 
|#(t)) = |tf*° (#(«))) and the fact that #|tf*°) = 



^2 m O\^m°)^m,n the auxiliary fields |$„) are related to 
the Born-Oppenheimer wave functions via Eq. A7 as 
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The auxiliary fields are inserted into the expre ssio n for 
the constraint forces acting on the atoms, Eq. A6, and 
after a few transformations the additional forces of the 
wave function cloud acting on the atoms are obtained: 
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The effective equation of motion given in Eq. for the 
nuclei is obtained by adding the corresponding constraint 
forces to the equations of motion for the nuclei. 
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